% KurmannOtrok_main.m

% Main program to replicate results in
%
%   Kurmann, A. and C. Otrok (2012). "News Shocks and the Slope of the Term Structure of Interest rates." American Economic Review (forthcoming)
%
% See readme.txt file for an explanation of the different programs.
%
% Written by Andre Kurmann, Federal Reserve Board. Last version: December 2012
% 
% Some small utilities used in the code are from J.P. Lesage's econometrics
% toolbox and are clearly acknowledged.
%
% Use of code for research purposes is permitted as long as proper reference to source is given. 
%
%------------------------------------------------------------------------------------------------------------------------------------------------

clear all
%close all
var_data=22; % variables included: 0=no investment, 1=no investment, no stock, no confidence , 11=full data, 10= no stock, 20=no stock+ no conf, 30=full+conf, 101=yfp,y,c,iv,pai, 110=tfp,y,c,iv,h,pai
cd_bands=1; %confiedence bands for fig 401
fig401=0; %fig400 is on if "1"
point_bays=0; % 1 if use point estimate for bayesisn

bps=0;
%-------------------------------------------------------------
% Step 0: parameter settings
%-------------------------------------------------------------
nlags = 4;                    % number of lags in VAR ; originally =4  
nconst = 1;                   % 1 = constant for , 0 otherwise :Even with no constant term estimation, 
                              % ** CHOOSE nconst=1 since we shut down a constant term in "three_var_final_hiro", "olsvarc_hiro" and related files 
                              % ** by forcing the estimate of a constant term to be zero (V matrics to be zero in "olsvarc_hiro").  
                             
prior=1;                      % for prior,  1 =  standard prior
                              %             inf = flat prior
ndraws=3000;                  % number of draws when computing error bands
bound=0.165;                   % 0.165 lower bound of confidence interval
                              % (e.g. for bound=0.16, CI is 16% - 84%)

nimp = 20;                    % horizon of VD and IRFs to be shown                        

KLbar=0;                      % lower horizon bound for max FEV share objective
KUbar=40;                     % higher horizon bound for max FEV share objective        
                                       
est_method = 2; % 0 : OLS point estimation
                % 1 : Bayesian estimation with Minnesota prior 
                % 2 : OLS point estimation and classical bootstrap             
                
shock_extract = 2; % 1: Uhlig's max FEV share shock
                   % 2: TFP news shock
                   % 3: Current TFP shock
                   
slope = 0;  % 1: maximize FEV share of slope if Uhlig method is used
            % 0: maximize FEV share of level if Uhlig method is used
              
use_yields = 0; % 0: use no yield curve variables
                % 1: use spread between long and 3-month rate and long-rate
                % 2: use yield spread b/w long and ffr -> long rate is not
                % in VAR but computed from spread and ffr thereafter
                % 3: use Diebold-Li slope and level factor

%-----------------------------------------------------------------------
% Step 1: load and transform data
%-----------------------------------------------------------------------
%esty1=1960;     % First Year for estimation     (max sample: 1959:2 - 2005:2 except if 10-year yield is used in which case sample starts in 1971:4)
%estq1=1;        % First quarter for estimation  

esty1=1957;     % First Year for estimation     
estq1=2;        % First quarter for estimation  
esty2=2004;     % Last Year for estimation      
%esty2=1999;     % Last Year for estimation      
estq2=4;        % Last quarter for estimation 

%three_var_final_hiro_percent; %percentage adjusted data
three_var_final_hiro_restat;
[y, x, vars, lev]=dataset_quarterly_sw_restat(esty1,estq1,esty2,estq2,nlags,shock_extract,use_yields,slope); 
%[y, x, vars, lev]=dataset_quarterly_sw_percent(esty1,estq1,esty2,estq2,nlags,shock_extract,use_yields,slope); 

%y is T x nvars matrix of rhs variables
%x is T x nvars*nlags of lhs variables
[T,nvars]=size(y);

%---------------------------------------------------------------------
% Step 3: Estimate VAR, extract shocks, FEV shares explained, and IRFs 
%---------------------------------------------------------------------
if nconst==1;
   const = ones(T,1); 
   x = [x const];
   ncoeffs = nvars*nlags+1;
else
   ncoeffs = nvars*nlags; 
end;

%ncoeffs + 1 x nvars matrix of regression coefficients
%b=x\y;       %OLS coefficients: ncoeffs x nvars (** more efficient way of computing least squares; i.e. b=inv(x'*x)*x'*y;
%b=inv(x'*x)*x'*y;

bb=x\y;%hiro 
%bb_const=bb(:,nvars);
%three_var_final_hiro;

b=beta_ols;

res = y - x*b;          %T x nvar matrix of residuals
vmat = (1/T)*(res'*res);
stds=sqrt(diag(vmat));

%OLS point estimates
if est_method==0; 
    if (shock_extract == 1)
        [data,v]=uhlig(b,res,vmat,nvars,nlags,KLbar,KUbar,nimp,slope,use_yields);
    elseif (shock_extract == 2)
        [data,v]=tfpnews(b,res,vmat,nvars,nlags,KLbar,KUbar,nimp,slope,use_yields);
        [data_k,v]=tfpnews(b,res,vmat,nvars,nlags,KLbar,KUbar,nimp,slope,use_yields);
                 datamed=data_k; 
        
    elseif (shock_extract == 3)
        [data,v]=current_tfp(b,res,vmat,nvars,nlags,nimp,slope,use_yields);
    end
    
%Bayesian estimation with Minnesota prior
elseif est_method==1;
     
    [hm,bm] = minneprc(y,x,nlags,1,nconst,lev,prior); %Minnesota prior
    xxx = inv(kron(eye(nvars),x'*x) + diagrv(eye(size(hm,1)),hm));
    bb = xxx * (vec(x'*y) + bm);    %stacked posterior means of regression coeffs: vec(bpost) = inv[inv(H) + kron(I,X'X)] * [vec[X'Y] + inv(H)*bprior]
    b=reshape(bb,ncoeffs,nvars);    %posterior coefficient matrix 
    res = y - x*b;                  %T x nvar matrix of residuals
    vmat = (1/T)*(res'*res);
    
    if ndraws==0;
        if (shock_extract == 1)
            [data,v]=uhlig(b,res,vmat,nvars,nlags,KLbar,KUbar,nimp,slope,use_yields);
        elseif (shock_extract == 2)    
            [data,v]=tfpnews(b,res,vmat,nvars,nlags,KLbar,KUbar,nimp,slope,use_yields);
             if point_bays==1
                 [data_k,v]=tfpnews(b,res,vmat,nvars,nlags,KLbar,KUbar,nimp,slope,use_yields);
                 datamed=data_k;
             end
        elseif (shock_extract == 3)
            [data,v]=current_tfp(b,res,vmat,nvars,nlags,nimp,slope,use_yields);
        end

    else
        sxx=chol(xxx)';
        sinv=chol(inv(vmat));
        datadraws=[];
        vdraws=[];
        randn('state',0);
        randmatrix=randn(nvars,T,ndraws);
        randvec=randn(nvars*ncoeffs,ndraws);
        for j=1:ndraws;
            [bbdraw,vmatdraw] = niw(b,sxx,sinv,T,randmatrix(:,:,j),randvec(:,j)); %drawing from posterior coefficient matrix
            bdraw=reshape(bbdraw,ncoeffs,nvars);
            res = y - x*bdraw;
            %compute results for each draw
            if (shock_extract == 1)
                [output,v]=uhlig(bdraw,res,vmatdraw,nvars,nlags,KLbar,KUbar,nimp,slope,use_yields);
            elseif (shock_extract == 2)              
                [output,v]=tfpnews(bdraw,res,vmatdraw,nvars,nlags,KLbar,KUbar,nimp,slope,use_yields);
            elseif (shock_extract == 3)
                [output,v]=current_tfp(bdraw,res,vmatdraw,nvars,nlags,nimp,slope,use_yields);
            end 
            datadraws(:,:,j) = output;
            vdraws(:,j) = v;
        end
        %compute median and confidence intervals   
        datadraws=sort(datadraws,3);
        low=round(ndraws*bound); high=round(ndraws*(1-bound));
        datalow=datadraws(:,:,low);
        %datamed=median(datadraws,3);
        datamed=mean(datadraws,3);
        datahigh=datadraws(:,:,high);
        data(:,:,1)=datamed; data(:,:,2)=datalow; data(:,:,3)=datahigh;
        vmed = median(vdraws,2);   %median values of shock
    end

%OLS and classical bootstrap for CIs       
elseif est_method==2;
    datadraws=[];
    naccept=0;    %counter for number of accepted sdraws
    for j=1:ndraws;
        %bnew=bootb(b,x,res,nlags,nconst);
        bnew=beta_boots(:,:,j);
        
        %bnew_=bootb_hiro(b,x,res,nlags,nconst);%hiro
        %BIAS_=[BIAS V]';%hiro
        %BIAS_=BIAS_(:,1:q);%hiro
        %bnew=bnew_+ BIAS_;%hiro
        
        res = y - x*bnew;          
        vmat = (1/T)*(res'*res);
        %compute results for each draw       
        if (shock_extract == 1)
            [output,v]=uhlig(bnew,res,vmat,nvars,nlags,KLbar,KUbar,nimp,slope,use_yields);
        elseif (shock_extract == 2)
            [output,v]=tfpnews(bnew,res,vmat,nvars,nlags,KLbar,KUbar,nimp,slope,use_yields);
        elseif (shock_extract == 3)
            [output,v]=current_tfp(bnew,res,vmat,nvars,nlags,nimp,slope,use_yields);
        end
            datadraws(:,:,j) = output;
            vdraws(:,j) = v;
        end
        %compute median and confidence intervals   
        datadraws=sort(datadraws,3);
        low=round(ndraws*bound); high=round(ndraws*(1-bound));
        datalow=datadraws(:,:,low);
        datamed=median(datadraws,3);
        %datamed=mean(datadraws,3);
        datahigh=datadraws(:,:,high);
        data(:,:,1)=datamed; data(:,:,2)=datalow; data(:,:,3)=datahigh;
        vmed = median(vdraws,2);   %median values of shock
end

%----------------------------------------------------------------------
% Step 3: Plots
%----------------------------------------------------------------------

%FEV shares explained by shock
%plotFEV(data,vars,slope,use_yields,shock_extract);
   
%IRFs to shock
plotIRF(data,vars,slope,use_yields,shock_extract)

% decomposition of yield curve response into EH part and term premia part
if use_yields == 1 | use_yields == 2
    plotEH(data,vars,slope,use_yields,shock_extract)
end


if fig401==1

figure(401)
%figure('NumberTitle','off','Name','IRFs of Monte Carlo VAR:  mean')

%figure('NumberTitle','off','Name','Monte carlo IRFs: VAR(median) & theoretical(median)')

if var_data==20
subplot(3,3,1)
plot(datalow(:,8),':'); hold on;  plot(datahigh(:,8),':'); hold on;  plot(datamed(:,8),'-'); 
title('tfp');
hold on

subplot(3,3,2)
plot(datalow(:,9),':'); hold on;  plot(datahigh(:,9),':'); hold on;  plot(datamed(:,9),'-'); 
title('y');
hold on

subplot(3,3,3)
plot(datalow(:,10),':'); hold on;  plot(datahigh(:,10),':'); hold on;  plot(datamed(:,10),'-'); 
title('c');
hold on

subplot(3,3,4)
plot(datalow(:,11),':'); hold on;  plot(datahigh(:,11),':'); hold on;  plot(datamed(:,11),'-'); 
title('iv');
hold on

subplot(3,3,5)
plot(datalow(:,12),':'); hold on;  plot(datahigh(:,12),':'); hold on;  plot(datamed(:,12),'-'); 
title('h');
hold on


subplot(3,3,6)
plot(datalow(:,13),':'); hold on;  plot(datahigh(:,13),':'); hold on;  plot(datamed(:,13),'-'); 
title('i');
hold on

subplot(3,3,7)
plot(datalow(:,14),':'); hold on;  plot(datahigh(:,14),':'); hold on;  plot(datamed(:,14),'-'); 
title('pai');
hold on

%subplot(3,3,8)
%plot(datalow(:,17),':'); hold on;  plot(datahigh(:,17),':'); hold on;  plot(datamed(:,17),'-'); 
%title('stock');
%hold on

end



if var_data==21
subplot(3,3,1)
plot(datalow(:,8),':'); hold on;  plot(datahigh(:,8),':'); hold on;  plot(datamed(:,8),'-'); 
title('tfp');
hold on

subplot(3,3,2)
plot(datalow(:,9),':'); hold on;  plot(datahigh(:,9),':'); hold on;  plot(datamed(:,9),'-'); 
title('y');
hold on

subplot(3,3,3)
plot(datalow(:,10),':'); hold on;  plot(datahigh(:,10),':'); hold on;  plot(datamed(:,10),'-'); 
title('c');
hold on

%subplot(3,3,4)
%plot(datalow(:,12),':'); hold on;  plot(datahigh(:,12),':'); hold on;  plot(datamed(:,12),'-'); 
%title('iv');
%hold on

subplot(3,3,5)
plot(datalow(:,11),':'); hold on;  plot(datahigh(:,11),':'); hold on;  plot(datamed(:,11),'-'); 
title('h');
hold on


subplot(3,3,6)
plot(datalow(:,12),':'); hold on;  plot(datahigh(:,12),':'); hold on;  plot(datamed(:,12),'-'); 
title('i');
hold on

subplot(3,3,7)
plot(datalow(:,13),':'); hold on;  plot(datahigh(:,13),':'); hold on;  plot(datamed(:,13),'-'); 
title('pai');
hold on

%subplot(3,3,8)
%plot(datalow(:,17),':'); hold on;  plot(datahigh(:,17),':'); hold on;  plot(datamed(:,17),'-'); 
%title('stock');
%hold on

subplot(3,3,9)
plot(datalow(:,14),':'); hold on;  plot(datahigh(:,14),':'); hold on;  plot(datamed(:,14),'-'); 
title('confidence');
hold on
end



if var_data==22
subplot(3,3,1)
plot(datalow(:,7),':'); hold on;  plot(datahigh(:,7),':'); hold on;  plot(datamed(:,7),'-'); 
title('tfp');
hold on

subplot(3,3,2)
plot(datalow(:,8),':'); hold on;  plot(datahigh(:,8),':'); hold on;  plot(datamed(:,8),'-'); 
title('y');
hold on

subplot(3,3,3)
plot(datalow(:,9),':'); hold on;  plot(datahigh(:,9),':'); hold on;  plot(datamed(:,9),'-'); 
title('c');
hold on

%subplot(3,3,4)
%plot(datalow(:,12),':'); hold on;  plot(datahigh(:,12),':'); hold on;  plot(datamed(:,12),'-'); 
%title('iv');
%hold on

subplot(3,3,5)
plot(datalow(:,10),':'); hold on;  plot(datahigh(:,10),':'); hold on;  plot(datamed(:,10),'-'); 
title('h');
hold on


subplot(3,3,6)
plot(datalow(:,11),':'); hold on;  plot(datahigh(:,11),':'); hold on;  plot(datamed(:,11),'-'); 
title('i');
hold on

subplot(3,3,7)
plot(datalow(:,12),':'); hold on;  plot(datahigh(:,12),':'); hold on;  plot(datamed(:,12),'-'); 
title('pai');
hold on

end

end


if cd_bands==1

figure(400)
%figure('NumberTitle','off','Name','IRFs of Monte Carlo VAR:  mean with percentiles')

%figure('NumberTitle','off','Name','Monte carlo IRFs: VAR(median) & theoretical(median)')


if var_data==20

subplot(3,3,1)
plot(datalow(:,8),':'); hold on;  plot(datahigh(:,8),':'); hold on;  plot(datamed(:,8),'-'); 
title('tfp');
hold on

subplot(3,3,2)
plot(datalow(:,9),':'); hold on;  plot(datahigh(:,9),':'); hold on;  plot(datamed(:,9),'-'); 
title('y');
hold on

subplot(3,3,3)
plot(datalow(:,10),':'); hold on;  plot(datahigh(:,10),':'); hold on;  plot(datamed(:,10),'-'); 
title('c');
hold on

subplot(3,3,4)
plot(datalow(:,11),':'); hold on;  plot(datahigh(:,11),':'); hold on;  plot(datamed(:,11),'-'); 
title('iv');
hold on

subplot(3,3,5)
plot(datalow(:,12),':'); hold on;  plot(datahigh(:,12),':'); hold on;  plot(datamed(:,12),'-'); 
title('h');
hold on


subplot(3,3,6)
plot(datalow(:,13),':'); hold on;  plot(datahigh(:,13),':'); hold on;  plot(datamed(:,13),'-'); 
title('i');
hold on

subplot(3,3,7)
plot(datalow(:,14),':'); hold on;  plot(datahigh(:,14),':'); hold on;  plot(datamed(:,14),'-'); 
title('pai');
hold on

end


if var_data==21
subplot(3,3,1)
plot(datalow(:,8),':'); hold on;  plot(datahigh(:,8),':'); hold on;  plot(datamed(:,8),'-'); 
title('tfp');
hold on

subplot(3,3,2)
plot(datalow(:,9),':'); hold on;  plot(datahigh(:,9),':'); hold on;  plot(datamed(:,9),'-'); 
title('y');
hold on

subplot(3,3,3)
plot(datalow(:,10),':'); hold on;  plot(datahigh(:,10),':'); hold on;  plot(datamed(:,10),'-'); 
title('c');
hold on

%subplot(3,3,4)
%plot(datalow(:,12),':'); hold on;  plot(datahigh(:,12),':'); hold on;  plot(datamed(:,12),'-'); 
%title('iv');
%hold on

subplot(3,3,5)
plot(datalow(:,11),':'); hold on;  plot(datahigh(:,11),':'); hold on;  plot(datamed(:,11),'-'); 
title('h');
hold on


subplot(3,3,6)
plot(datalow(:,12),':'); hold on;  plot(datahigh(:,12),':'); hold on;  plot(datamed(:,12),'-'); 
title('i');
hold on

subplot(3,3,7)
plot(datalow(:,13),':'); hold on;  plot(datahigh(:,13),':'); hold on;  plot(datamed(:,13),'-'); 
title('pai');
hold on

%subplot(3,3,8)
%plot(datalow(:,17),':'); hold on;  plot(datahigh(:,17),':'); hold on;  plot(datamed(:,17),'-'); 
%title('stock');
%hold on

subplot(3,3,9)
plot(datalow(:,14),':'); hold on;  plot(datahigh(:,14),':'); hold on;  plot(datamed(:,14),'-'); 
title('confidence');
hold on
end


if var_data==22
subplot(3,3,1)
plot(datalow(:,7),':'); hold on;  plot(datahigh(:,7),':'); hold on;  plot(datamed(:,7),'-'); 
title('tfp');
hold on

subplot(3,3,2)
plot(datalow(:,8),':'); hold on;  plot(datahigh(:,8),':'); hold on;  plot(datamed(:,8),'-'); 
title('y');
hold on

subplot(3,3,3)
plot(datalow(:,9),':'); hold on;  plot(datahigh(:,9),':'); hold on;  plot(datamed(:,9),'-'); 
title('c');
hold on

%subplot(3,3,4)
%plot(datalow(:,12),':'); hold on;  plot(datahigh(:,12),':'); hold on;  plot(datamed(:,12),'-'); 
%title('iv');
%hold on

subplot(3,3,5)
plot(datalow(:,10),':'); hold on;  plot(datahigh(:,10),':'); hold on;  plot(datamed(:,10),'-'); 
title('h');
hold on


subplot(3,3,6)
plot(datalow(:,11),':'); hold on;  plot(datahigh(:,11),':'); hold on;  plot(datamed(:,11),'-'); 
title('i');
hold on

subplot(3,3,7)
plot(datalow(:,12),':'); hold on;  plot(datahigh(:,12),':'); hold on;  plot(datamed(:,12),'-'); 
title('pai');
hold on

%subplot(3,3,8)
%plot(datalow(:,17),':'); hold on;  plot(datahigh(:,17),':'); hold on;  plot(datamed(:,17),'-'); 
%title('stock');
%hold on

end



end





if cd_bands==0

figure(400)
%figure('NumberTitle','off','Name','IRFs of Monte Carlo VAR:  mean with percentiles')

%figure('NumberTitle','off','Name','Monte carlo IRFs: VAR(median) & theoretical(median)')


if var_data==20
subplot(3,3,1)
plot(datamed(:,8),'-'); 
title('tfp');
hold on

subplot(3,3,2)
plot(datamed(:,9),'-'); 
title('y');
hold on

subplot(3,3,3)
 plot(datamed(:,10),'-'); 
title('c');
hold on

subplot(3,3,4)
plot(datamed(:,11),'-'); 
title('iv');
hold on

subplot(3,3,5)
plot(datamed(:,12),'-'); 
title('h');
hold on


subplot(3,3,6)
plot(datamed(:,13),'-'); 
title('i');
hold on

subplot(3,3,7)
plot(datamed(:,14),'-'); 
title('pai');
hold on
end


if var_data==21
subplot(3,3,1)
plot(datamed(:,8),'-'); 
title('tfp');
hold on

subplot(3,3,2)
plot(datamed(:,9),'-'); 
title('y');
hold on

subplot(3,3,3)
plot(datamed(:,10),'-'); 
title('c');
hold on

%subplot(3,3,4)
%plot(datalow(:,12),':'); hold on;  plot(datahigh(:,12),':'); hold on;  plot(datamed(:,12),'-'); 
%title('iv');
%hold on

subplot(3,3,5)
plot(datamed(:,11),'-'); 
title('h');
hold on


subplot(3,3,6)
plot(datamed(:,12),'-'); 
title('i');
hold on

subplot(3,3,7)
  plot(datamed(:,13),'-'); 
title('pai');
hold on

%subplot(3,3,8)
%plot(datalow(:,17),':'); hold on;  plot(datahigh(:,17),':'); hold on;  plot(datamed(:,17),'-'); 
%title('stock');
%hold on

subplot(3,3,9)
 plot(datamed(:,14),'-'); 
title('confidence');
hold on
end


if var_data==22
subplot(3,3,1)
 plot(datamed(:,7),'-'); 
title('tfp');
hold on

subplot(3,3,2)
 plot(datamed(:,8),'-'); 
title('y');
hold on

subplot(3,3,3)
 plot(datamed(:,9),'-'); 
title('c');
hold on

%subplot(3,3,4)
%plot(datalow(:,12),':'); hold on;  plot(datahigh(:,12),':'); hold on;  plot(datamed(:,12),'-'); 
%title('iv');
%hold on

subplot(3,3,5)
 plot(datamed(:,10),'-'); 
title('h');
hold on


subplot(3,3,6)
  plot(datamed(:,11),'-'); 
title('i');
hold on

subplot(3,3,7)
 plot(datamed(:,12),'-'); 
title('pai');
hold on

%subplot(3,3,8)
%plot(datalow(:,17),':'); hold on;  plot(datahigh(:,17),':'); hold on;  plot(datamed(:,17),'-'); 
%title('stock');
%hold on

end

end




